Molecular characterization of recombinant LSDV isolates from 2022 outbreak in Indonesia through phylogenetic networks and whole-genome SNP-based analysis

Lumpy skin disease (LSD) is a transboundary viral disease of cattle and water buffaloes caused by the LSD virus, leading to high morbidity, low mortality, and a significant economic impact. Initially endemic to Africa only, LSD has spread to the Middle East, Europe, and Asia in the past decade. The most effective control strategy for LSD is the vaccination of cattle with live-attenuated LSDV vaccines. Consequently, the emergence of two groups of LSDV strains in Asian countries, one closely related to the ancient Kenyan LSDV isolates and the second made of recombinant viruses with a backbone of Neethling-vaccine and field isolates, emphasized the need for constant molecular surveillance. This current study investigated the first outbreak of LSD in Indonesia in 2022. Molecular characterization of the isolate circulating in the country based on selected LSDV-marker genes: RPO30, GPCR, EEV glycoprotein gene, and B22R, as well as whole genome analysis using several analytical tools, indicated the Indonesia LSDV isolate as a recombinant of LSDV_Neethling_vaccine_LW_1959 and LSDV_NI-2490. The analysis clustered the Indonesia_LSDV with the previously reported LSDV recombinants circulating in East and Southeast Asia, but different from the recombinant viruses in Russia and the field isolates in South-Asian countries. Additionally, this study has demonstrated alternative accurate ways of LSDV whole genome analysis and clustering of isolates, including the recombinants, instead of whole-genome phylogenetic tree analysis. These data will strengthen our understanding of the pathogens’ origin, the extent of their spread, and determination of suitable control measures required. Supplementary Information The online version contains supplementary material available at 10.1186/s12864-024-10169-6.


Introduction
Lumpy skin disease virus (LSDV) is a member of the Capripoxvirus (CaPV) genus that has a linear doublestranded DNA genome [1].The virus causes lumpy skin disease (LSD) in cattle and water buffaloes [2].The CaPV genus comprises two other species, the sheeppox virus (SPPV) and goatpox virus (GTPV), which are antigenically and genetically closely related to LSDV [1,3].LSD infections are associated with high morbidity (up to 100%) and low mortality (< 10%), with a significant economic impact due to losses related to animal production, restrictions to trade, and treatment/vaccination costs, and therefore a significant World Organisation for Animal Health (WOAH) notifiable disease [2,4].Clinically sick animals present symptoms such as fever, enlarged lymph nodes, nasal discharge, and nodular skin lesions [5,6].
LSDV is transmitted by blood-feeding vectors such as stable flies, mosquitoes, and ticks [7][8][9].This may have facilitated its spread from the endemic region in Africa to the Middle East, Europe, and South Asia in the past ten years [10].Other modes of transmission, such as direct contact, seminal, and intrauterine transmission, have also been demonstrated experimentally [11][12][13].LSD control measures include cattle movement/trade restriction, vector control, disinfection, or culling of the animals showing clinical signs [2,14], with the most effective control strategy being vaccination using live attenuated vaccines, such as those derived from the Neethling strain [15][16][17].
Understanding the molecular and genetic epidemiology of CaPVs is a prerequisite for better management and control of LSD, given the antigenic similarity of CaPVs that challenges their distinction by serological tests.Among several molecular-based techniques for genetically distinguishing CaPVs include, the sequence analysis of genes encoding the RNA polymerase 30 kDa subunit (RPO30), the G-protein-coupled receptor (GPCR), the CaPV homolog of the variola virus B22R, and the EEV glycoprotein that possess specific signatures that can differentiate the CaPVs [18][19][20][21].These genes could also help distinguish live attenuated vaccines from field isolates and the recombinant LSDV strains.However, a whole genome-based approach is preferable as it enables proper recombination analysis which has been done for the recently reported recombinants [22,23].Nevertheless, phylogenetic tree construction using recombinants' whole genome is not suitable/accurate since the analysis process considers mainly the mutations and speciation events and does not account for recombination between the sequences or ancestor sequences [24].Therefore, alternative tools such as restricting the analysis to parts of the virus genome, or phylogenetic networks and single nucleotide polymorphism (SNP) genotyping, which consider events like recombination/hybridization have been shown as more appropriate for whole genome analysis of recombinants [25][26][27].
LSD was first reported in Asia in 2019 in India, Bangladesh, and China and subsequently spread to other Asian countries such as Nepal, Myanmar, Thailand, and Vietnam [28].Phylogenetic analysis of the Asian LSDV isolates has demonstrated that the isolate circulating in South-Asian countries resembles the ancient LSDV isolate from Kenya [29][30][31], while the isolates in Southeast and Eastern Asia are recombinant viruses with Neethling vaccine strain and field isolates backbones [23,[32][33][34][35].In Indonesia, the initial suspected clinical cases of LSD were reported in November 2020 cattle and buffaloes from Kampar and buffaloes from Kuansing districts, Riau Province.In 2021, more suspicions were reported in cattle from Padang Pariaman district in North Sumatera (January) and in a buffalo from Muara Bungo district (August).However, none of the suspected LSD cases from 2020 to 2021 could be confirmed by either serological (Capripox ELISA) or molecular tests (Polymerase chain reaction (PCR).In February 2022, the disease spread to cattle from Indragiri Hulu district, where samples were collected and confirmed by PCR at the Research Centre for Veterinary Science (RCVS), National Research and Innovation Agency, enabling the notification to WOAH on February 2022 [36].This study aims to provide clinical evidence of the first outbreak of LSD in Indonesia in 2022 and the molecular characterization of the LSDV isolate circulating in the country based on selected LSDV-marker genes and whole genome analysis.Due to the recent emergence of LSDV recombinants that make it challenging (less accurate) to cluster LSDV isolates using whole genome phylogenetic trees, the study will also explore some alternative approaches based on existing tools for better accuracy in LSDV whole genome-based clustering.

LSD outbreak investigation and clinical signs
The suspected LSD infected animals presented nodular lesions on different body parts, nasal discharge, and fever (Fig. 1).The morbidity and mortality rates recorded in the Indragiri Hulu district by February 2022 were 28.04% (30/107) and 1.87% (2/107), respectively.

Molecular detection and diagnosis of LSDV
CaPV genome was detected in ten (83.33%) out of the twelve samples collected from 2022 suspected cattle cases in Indragiri Hulu district, Indonesia, by RT-PCR using the Bowden et al. protocol [6].Further analysis of the positive samples, using HRM assay, confirmed that the suspected cases were infected with LSDV (Additional file 1: Table s1).

Targeted gene sequence and phylogenetic analysis
For five samples, the targeted fragments were successfully amplified and sequenced, and the assembled contigs of 606, 1146, 327, and 832 bp for the RPO30, GPCR, EEV, and B22R genes, respectively, were deposited in the Gen-Bank with accession numbers OR246967 to OR246986.
The multiple sequence alignments showed that all five Indonesia LSDV samples were 100% identical based on RPO30, GPCR, EEV, and B22R gene sequences.Phylogenetic analysis of the complete RPO30 gene sequences showed that Indonesia sequences were grouped in LSDV Cluster I together with LSDV Neethling-derived vaccines and other LSDV isolates from East and Southeast Asia, such as China, Taiwan, Vietnam, and Thailand (Fig. 2A).However, the GPCR gene sequences grouped the Indonesia sequences in LSDV Cluster II with the common LSDV field isolates, the historical LSDVs Fig. 2 Neighbour-joining tree based on the complete (A) RPO30 and (B) GPCR gene sequences of CaPVs with LSDVs from Indonesia (in red) visualized on iTOL.The evolutionary distances were computed using the Maximum Composite Likelihood method with 1000 bootstrap replicates on MEGA X Fig. 1 Clinical signs of LSD suspected cases in Indonesia; nodules on cattle from Indragiri Hulu district from Kenya, and LSDV KSGP-0240 derived vaccine, yet with the other East and Southeast Asia isolates (Fig. 2B).
The multiple sequence alignment of the partial EEV glycoprotein gene showed that Indonesian LSDV isolates are identical in this region of the genome to LSDV isolates from China, Vietnam, Taiwan, and LSDV Neethling-derived vaccines that have a 27-nucleotide deletion (175-201) (Fig. 3A).The B22R sequence alignment showed that part of the Indonesian sequence was identical to the historical LSDVs from Kenya (353 bp) and the other to the LSDV Neethling-derived vaccine (479 bp).However, the B22R of the Indonesian LSDVs differed from LSDV Neethling and LSDV KSGP-0240 derived vaccines due to one nucleotide InDel at positions 102 and 745, respectively (Fig. 3B).
Similar recombination events are also present in other LSDV isolates from China, Thailand, Hong Kong, Vietnam, Taiwan, and Tomsk, Russia (Table 1).Some of the ORFs affected by the recombination events include the G-protein coupled receptor gene (LSDV011), the DNA polymerase gene (LSDV039), the DNA ligase-like protein gene (LSDV133), and the variola virus B22R-like protein gene (LSDV134).These four ORFs are among the nine proteins annotated as only 99.5-99.9%homologous to both LSDV_NI-2490 and LSDV_Neethling_vac-cine_LW_1959 reference genomes (Additional file 1: Table s3).The most notable ORF was LSDV087, whose sequence analysis showed a nucleotide insertion at position 599 in all Neethling vaccine strains which caused a frameshift and truncated the ORF to 603nt instead of 763nts.No recombination was observed within this ORF, and all known LSDV recombinants were 100% identical to the field isolates in this region of the genome.Further, in addition to the 8 SNPs observed between the LSDV NI-like field isolates and Neethling-like field isolates, the nucleotide alignment of ORF087 showed 14 SNPs between the LSDV NI-like field isolates and Neethling derived vaccines; and 7 SNPs between Neethling-like field isolates and Neethling derived vaccines (Additional file 2: Figure s1).
The phylogenetic analysis of sequence fragments located on the 5' and 3' ends of breakpoint events supported the recombination events.For example, analysis of the sequences located on the 5' end of the breakpoint event 10, clustered the Indonesian LSDVs with LSDV_ NI-2490 (Cluster II), while the sequences on the 3' end of the breakpoint clustered them with LSDV_Neeth-ling_vaccine_LW_1959 (Cluster I) (Additional file 3: Figure s2).The recombination of Indonesian LSDVs was also confirmed using SimPlot, followed by phylogenetic analysis of the sequences located on the 5' and 3' ends of the breakpoints (Additional file 4: Figure s3).Similarly, the phylogenetic network analysis of the aligned LSDV genomes using PopART software indicated that the Indonesia LSDVs are likely recombinants of NI-like and Neethling-like LSDVs placed on the opposite ends of the network.The network confirmed that the LSDVs from Indonesia cluster together with other isolates from East and Southeast Asia (recomb_3-like), which have been demonstrated as LSDV recombinants and are different from the known LSDV recombinants from Russia (recomb_1-like, recomb_2-like, and recomb_4-like) (Fig. 5).

Discussion
This study has presented typical clinical signs of LSD including skin nodules, nasal discharge, and high fever in suspected cattle cases in 2022 outbreak in Indragiri Hulu district, Indonesia, with morbidity and mortality rates of  28.04% and 1.87%, respectively.It has been reported that the severity of LSD varies based on the immune status of the animal, the breed, the age, and the virus strain, leading to an average mortality rate of < 10% and variable morbidity rate from 5 to 45% [2,4].However, the morbidity and mortality rates recorded in the LSD outbreak in Indonesia were consistent with other LSD outbreaks in Asian countries [28,37,38].
The initial molecular characterization of the Indonesian LSDV isolates based on targeted CaPV-marker genes indicated possible virus recombination [18][19][20][21].First, the nucleotide sequence analysis of the four targeted genes showed that the Indonesian LSDV isolates were 100% identical to the recently demonstrated LSDV recombinants from East and Southeast Asia: China, Thailand, Taiwan, Vietnam, and HongKong [23,[32][33][34][35]. Second, on one hand, the complete RPO30 gene and the partial EEV gene sequence analysis showed that Indonesian LSDV sequences have 99.67% and 100% sequence similarity, respectively, to the LSDV Neethling-derived vaccines.On the other hand, the GPCR gene sequence similarity was highest (99.13%) to the common LSDV field strains, while the B22R sequence was 99.40% similar to the Neethling vaccine strain and 99.28% to the NI-2490 LSDV strain.
Further analysis of the Indonesian LSDV whole genome sequences confirmed its closest homology to the recombinant LSDVs from East and Southeast Asia and that it was 99.56% similar to Neethling_LSDV_vaccine and 99.25% to the LSDV_NI-2490 strain.Possible recombination of the Indonesian LSDV genomes was assessed and verified using different tools (RDP4, SimPlot, DAPCA, and PopART network) that have previously been applied in recombination analysis [39,40].The PopART network and the PCA scatter plot showed seven clusters of LSDVs and grouped Indonesia isolates with other reported East and Southeast Asia recombinant LSDVs [41].The RDP4 and SimPlot indicated that Indonesia LSDVs are most likely recombinants of Neethling_LSDV_vaccine and LSDV_NI-2490 isolates with 11 recombination events.Notably, the recombinant events detected in this study have been reported in LSDV genomes from Thailand and China [32,34].Although the parental donors for all the reported LSDV recombinants are Neethling vaccine and Kenyan NI-2490-like strains, the Indonesia LSDV recombination was similar to other recombinants from East and Southeast Asia but different from those found in Russia that have 27 (LSDV_ Russia/Saratov/2017), 24 (LSDV_Udmurtiya_Russia_2019) and 25 (Tyumen_2019, OL542833_1) recombination events [42][43][44].The difference in recombination patterns of LSDVs was demonstrated in this study by the SNP heatmap and the PopART network analysis that showed the phylogenetic relatedness of the recombinants to their parental donors.
In addition to the previously reported genes employed in combination in this study (GPCR, RPO30, EEV, and B22R), the LSDV whole genome analysis suggested that LSDV087 could be an additional DIVA target through PCR and sequencing, that can complement the existing CaPV RT-PCR DIVA methods, due to the unique signatures or SNPs found at DNA level in LSDV Neethling derived vaccines compared to the field strains including the LSDV recombinants [18][19][20][21].In addition to the intergenic recombination (exchange of full-length gene sequence) observed in Indonesian LSDV genomes and other reported recombinant LSDVs, intragenic recombination (exchange of gene sequence fragments) was also detected in nine genes [45].Recombination in viruses is common in the field following a co-infection of viruses in a single cell of a host which is influenced by the selection pressure of a virus population and contributes to virus evolution through genetic shift [45,46].The intragenic recombinations like those detected in the Indonesia LSDV isolates are related to the regulation of gene expression, hence altering the pathogenic properties of the virus [45,46].Virus recombination in the field is likely following vaccination with live attenuated vaccine in a population where wild-type virus is already circulating or contamination in cell culture during vaccine production [47].No LSDV vaccination program was reported in Indonesia before the LSD outbreak; therefore, it is possible that the virus was introduced into the country from the neighboring countries.This is likely considering that the Indonesian LSDV isolate is nearly identical to the isolates circulating in the neighboring countries in East and Southeast Asia, suggesting a common origin of the virus in the region.A recent study has shown that the emergence of LSDV recombinants in Asia is most likely a result of a spillover from animals vaccinated with the Lumpivax vaccine that appears to be a mixture of Neethling-like LSDV vaccine strain, KSGP-0240 vaccine strain and several LSDV recombinants that are nearly identical to those circulating in Asia [41].For instance, the initial outbreak of the recombinant_3-like LSDVs was reported in Qapqal Xibe county, Xinjiang province in Northwest China, which is located ∼ 20 km from the border of Kazakhstan, where a mass vaccination campaign with Lumpivax vaccine was launched [48].Thus, it is likely that in addition to the poor-quality vaccine spillover, cattle movement, trade and insect vectors may have contributed to the spread of the LSDV recombinants into the neighboring countries in the region.Moreover, It would be interesting to further investigate whether the changes in proteins for viral virulence and host range, such as G-protein coupled receptor (LSDV011), IFN-alpha beta receptor glycoprotein (LSDV135), IL-10 (LSDV005) and IL-1 receptor (LSDV006), that were found in LSDV recombinants circulating in Indonesia and other Southeast Asian countries cause higher infectivity and virulence, thereby contributing to the virus's rapid spread as suggested in earlier studies with LSDV recombinants in Russia [49,50].

Conclusion
This study has reiterated the importance of applying several analytical tools for accurate LSDV whole genome analysis, especially with datasets including recombinants.Based on these tools, the study has shown that the Indonesia LSDV isolate is a recombinant of Neeth-ling_vaccine_LW_1959 and LSDV_NI-2490 strains and is closely related to the previously reported LSDV recombinants circulating in East and Southeast Asia, suggesting a common origin/source of the LSDV isolate causing the outbreak and spreading in the region.Most importantly, the analysis in this study has identified a unique ORF, LSDV087, as a potential DIVA target for LSDVs through sequencing.Therefore, this study demonstrates the significance of accurate molecular characterization of LSDV for disease surveillance and understanding of its origin, as well as vaccine quality control before animal vaccination.

Sample collection and processing
Six skin nodules and six nasal swab samples were collected from suspected cattle cases in Indragiri Hulu district in 2022 and sent to the RCVS for further investigation.DNA was extracted using the Qiagen AllPrep DNA/RNA extraction kit and eluted in 50 μl elution buffer.

PCR and differential diagnosis
A real-time PCR (RT-PCR), previously described by Bowden et al. [6], was used to detect the CaPV genome in all twelve collected samples.The PCR reactions were prepared in a 20 μl reaction volume containing 18 μl of the PCR master mix with primers, a probe, and two μl of template DNA.The positive samples were characterized further using the high-resolution melting (HRM) assay capable of differentiating CaPVs based on the melting temperature (Tm) of the PCR amplicon [51].

Amplification and sequencing of selected CaPV genes
Four selected CaPV-marker genes (RPO30, GPCR, EEV, and B22R) that can differentiate LSDV vaccine strains from field isolates were amplified and sequenced for the CaPV-positive samples following the previously described PCR conditions [52].The PCR products were analyzed on a 2% agarose gel in 1X TAE buffer and visualized on a Gel Documentation System (Bio-Rad).The PCR products were purified and sent for Sanger sequencing from both directions using forward and reverse primers at LGC Genomics (Germany).

Whole genome sequencing of LSDV using Ion S5 technology
Approximately 100 ng of DNA was enzymatically fragmented into 200 bp lengths using Ion shear Plus reagents.The fragmented DNA was ligated to adapters and barcodes using the Ion Xpress™ Plus Fragment Library Kit and the Ion Xpress barcode adapters (Thermo Fisher Scientific, USA) following the manufacturer's instructions.Size selection was performed using Pippin Prep (Sage Science, Inc., USA).The libraries were further amplified for eight cycles using Platinum™ PCR SuperMix high fidelity and library amplification primer mix provided with the Ion Plus Fragment Library Kit.The amplified barcoded libraries at a concentration of 100 pM were pooled in equal volumes and loaded onto the Ion Chef™ Instrument (Thermo Fisher Scientific, USA) for automated template preparation and chip loading using Ion 540™ Kit-Chef (Thermo Fisher Scientific, USA).Using the Ion Chef™ Instrument, the pooled libraries were clonally amplified on the Ion spheres (ISPs) by emulsion PCR, followed by automated loading template-enriched ISPs onto an Ion 540 chip.Sequencing was performed on an Ion S5™ Next generation sequencing system (Thermo Fisher Scientific, USA) with 500 flows to generate 200 bp reads.

Sequence and phylogenetic analysis of CaPV-marker genes
The obtained sequences were assembled using Vector NTI software (Invitrogen) v11.5.For each targeted gene, the sequences were analyzed with those of other LSDVs, SPPVs, and GTPVs retrieved from GenBank.Multiple sequence alignments were performed on MEGA X using the Muscle algorithm and the codon option.Neighborjoining trees for the complete RPO30 and GPCR genes were also built on MEGA X, and the evolutionary distances were computed using the Maximum Composite Likelihood method with 1000 bootstrap replicates [53].The phylogenetic trees with the associated metadata were visualized on the Interactive Tree of Life (ITOL) tool [54].The multiple sequence alignments, analysis, and visualization of the partial EEV glycoprotein and B22R genes were conducted using BioEdit v7.2.5.
Given that the Indonesia LSDVs were grouped in LSDV Cluster I or II based on different LSDV-marker genes analysis, the open reading frames (ORFs) of the assembled Indonesia LSDV genomes were predicted with GATU using reference genomes from LSDV Cluster I; LSDV_Neethling_vaccine_LW_1959 (AF409138.1)and Cluster II; LSDV_NI-2490 (AF325528.1)[55].The whole genome sequences of Indonesia LSDV isolates, LSDV_ Indonesia_2022_S1, and LSDV_Indonesia_2022_S4, were deposited in the GenBank database under the accession numbers OR232413 and OR232414, respectively.

LSDV whole genome SNP and recombination analysis
The assembled Indonesian LSDV genomes were aligned with 46 complete LSDV genome sequences retrieved from GenBank using MAFFT.The SNPs from the aligned LSDV dataset relative to the LSDV_NI-2490 (AF325528.1)genome were analyzed using the adegenet package in R. First, the alignment file was converted into a genlight object using the function fasta2genlight and the SNPs were extracted [56].The SNPs distribution and density of the LSDV genotypes were visualized by heatmap.The discriminant analysis of Principal Component Analysis (DAPCA) of the extracted SNPs was also performed to identify the different clusters of the LSDV genotypes, which were plotted as a PCA scatterplot and neighbor-joining tree.Further SNP-based analysis was performed on all LSDV genomes clustering together with the Indonesia LSDV isolates.Three different tools were applied to assess evidence of recombination in the Indonesia LSDV isolate.The potential recombination events of the Indonesia LSDV isolates were first evaluated using the RDP4 software, which includes seven methods (RDP, GENECONV, BootScan, MaxChi, Chimera, SiScan, and 3Seq) [57].The recombination was then analyzed using the SimPlot 3.5.1 software with a window size of 1000 bp and a step of 20 bp.Given that recombinant genome sequences are not suitable for whole genome phylogenetic analysis [25,27], several phylogenetic trees were constructed based on the fragment sequences before and after the breakpoints predicted by RDP4 and Sim-Plot software on MEGA X, using the neighbor-joining method with 1000 bootstraps and visualized on ITOL.Lastly, the recombination was explored by conducting a phylogenetic network analysis of the aligned LSDV genomes on the PopART program, using the medianjoining method with epsilon set to zero [40].

Fig. 3
Fig. 3 Multiple sequence alignment of (A) the partial sequences of the EEV glycoprotein gene showing the 27-nucleotide deletion highlighted in the block found in Indonesia isolates and (B) partial sequences of the B22R gene showing the nucleotide insertion (in blocks) in LSDV_Neethling and LSDV_KSGPO-240 vaccines that are absent in Indonesia isolates (in red).The dots indicate the identical nucleotides in the alignment

Fig. 4
Fig. 4 Analysis of SNPs extracted from LSDV genome alignments.(A) SNP distribution heatmap, (B) PCA-scatter plot, and neighbor-joining tree based on the PCA scores of the SNP data of all LSDV complete genomes, (C) SNP distribution heatmap of recombinant_3-like LSDV genomes

Fig. 5
Fig. 5 Median-joining network inferred from LSDV whole genome sequences using the PopART program.The network shows the Indonesia LSDVs (in red) in cluster Recomb_3-like, between the NI-like and Neethling-like LSDV clusters.The number of mutations between each genome is labeled

Table 1
Predicted potential recombination events using different detection methods available in the RDP4 program *= The actual breakpoint position is undetermined (it was most likely overprinted by a subsequent recombination event) NS = Not significant